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Monte Carlo simulation using a cluster algorithm is used to compute the scaling 
part of the free energy for a three dimensional 0(4) spin model. The results are 
relevant for analysis of lattice studies of high temperature QCD. 



Introduction 



The high temperature phase transition for QCD with two flavors of hght quarks is ex- 
pected to be driven by chiral symmetry restoration, with an order parameter having 0(4) 
symmetry in the continuum hmitjl], ^. Thus, near the transition we expect the scahng 
properties of a three dimensional 0(4) spin model. For quark mass or temperature not too 
close to the transition, the system would be expected to behave like mean field theory. Re- 
cently Kogut and Kocic have suggested that mean field behavior might describe the system 
arbitrarily close to the critical point Finally, with Kogut-Susskind quarks on a nonzero 
lattice spacing, the exact chiral symmetry is only 0(2), and it is possible that lattice sim- 
ulations are better described by 0(2) critical behavior. In addition to its intrinsic interest 
as an indicator of the physics of the transition, the form of the free energy near the critical 
point is important in extrapolating the QCD equation of state from the quark masses where 
lattice simulations are practical to the light quark masses of the real world 0. 

Assuming a second order transition0, we expect the singular behavior of thermodynamic 
observables near the transition to be universal, meaning that the symmetry group of the order 
parameter and the dimension of the system determine the critical exponents and the form 
of the singular part of the free energy, up to normalization of the scaling variables. (See, for 
example, [§]) The critical exponents for 0(4) and 0(2) are well knownP, 0, but the form 
of the free energy, or the "scaling function", is only poorly known. An epsilon expansion 
result is available ||12||, quoted in Ref. 0. Similarly, Monte Carlo calculations of critical 



exponents have been used to study the critical behavior of high temperature QCD|l3[, but 
to date the full power of the scaling ansatz, namely comparison with the universal scaling 
functions as well as critical exponents, has not been brought to bear. 

Here we use Monte Carlo simulation to compute an approximate scaling function for 
0(4), to be used in comparing to Monte Carlo simulations of QCD. 

For the 0(N) spin model we use the partition function: 



Z = J[ds\exp ijJ^Si ■ Sj + soi J 



where ij are nearest neighbor pairs on a (hyper) cubic lattice in d dimensions and soi is the 
zero component of Sj. Then the energy and magnetization are: 

In QCD, using the normalization where the plaquette (□) is three when all links are unity. 
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Figure 1: A length rescaling by a factor of b is accomplished by changing the couplings (t, h) 
at point A to {byH,by'^h) at point B by moving along the renormalization group trajectory 
(curved line). The trajectories may be labelled by their intersections with the unit circle, 
so specifying the free energy on the unit circle, together with the scaling ansatz which tells 
how the free energy changes along a trajectory, specifics the free energy everywhere. The 
discontinuity in the order parameter at t < and h — (heavy line) implies that the 
derivative of the free energy is discontinuous there. 



the analogous equations are 

1 dlog{Z) 



2Vnt 



J 1 glog(Z) 

^^^^ = v^,^^ 



Parameterizing the scaling functions 

From invariance under a length rescaling by a factor b, the critical part of the free energy 
should have the property: 

fs{t,h)^b-''fs{h^%by^h) . (4) 

Here t and h are the scaling variables, with the critical point at (t, h) = (0, 0), and yt and 
Uh are the corresponding critical exponents. Other exponents can be expressed in terms of 
yt and yh- t = {T — Tc)/Tq and h = H/Hq are conventionally normalized by requiring that 
M{t = Q,h) = h}^^ and M{t < 0, /i = 0) = {—t)^- The free energy also has a nonsingular 
part. 
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The scaling ansatz, Eq. ^ implies that the magnetization near the critical point is deter- 
mined by a universal scaling function, conventionally written as: 

^ = f{t/h'/^') = f{x) . (5) 

The normalization conditions on t and h then require that /(O) = 1 and f{x) — > (— x)^ as 
X — oo. 



In computing the energy and pressure of QCD, we require the plaquette, analogous to 
the energy in the spin model, extrapolated to zero quark mass, which is analogous to zero 
magnetic field. The magnetization (or tptp in QCD) is y ^ ^°q}^^ = Yh^ ^ ^°ah^^ while the en- 
ergy (or plaquette in QCD) is y^^^^ = v% ^'°at^^ - Since the energy and magnetization 
are derivatives of the free energy with respect to t and h respectively, information about 
one quantity constrains the other. In particular, the behavior of ipip in QCD can help in 
extrapolating the plaquette to zero or small quark mass. (It is important to enforce consis- 
tency of the plaquette and ipip in computing the equation of state. Using a plaquette and 
ipijj which are not derivatives of the same free energy could lead to inconsistent thermody- 
namics.) Because our analysis requires that we handle the energy and magnetization on the 
same footing, in addition to Eq. |^, we develop a formulation of the scaling ansatz which 
treats the magnetization and energy equally. 

The scaling ansatz, Eq. |^, tells us that if we specify the singular free energy on any circle 
in the t,h plane, we have specified it for all t,h. (See Fig. |ip In particular, if g{6) is the 
scaling free energy on the unit circle in the t, h plane, then the rescaling factor b which takes 
t, h to the unit circle is determined by 

(6^*t)^ + {by^hf = 1 (6) 

For yt and i/h positive this clearly has a unique solution for 6 > given t and h. Although 
in general this cannot be solved analytically for b, it is straightforward to do differentiations 
implicitly and solve the equation numerically. Then the singular free energy (actually minus 
one times temperature times free energy per volume) is 

^\ogiZ,{t,h)) = b{t,h)-''g{9it,h)) (7) 

where 6{t, h) = atan2(6^'"/i, 6^*t), and g{6) is a universal function. 

Equivalently, the relation between t, h and b, 9 can be expressed as: 

b^^h = sm{e) 

= cos{e) . (8) 
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After some differentiations, the magnetization and energy can be expressed in terms of 



where 



db. -hb'^y^ - sm(e)by'^+^ 



db. -tb'^y -cosieW'^^ 



dt yttWy^-^ + yhhWy^-^ yt cos^{e) + yh sm\9) 



ah" dh"" "•' y,ms^(e)+y^sia'(e) 

^1 - -hb--y^ + ht'^ (y, - y,) 6---- = Slf^:^'}., 



dt''' yt cos^ (6) +yhSin' (6) 



(11) 



Physical insight into the form of g{6) comes from considering special cases. 
First, for i > and h small: 

Ho 

E^9(0)^d/yt-i (12) 
Toyt ^ ' 

Since M must vanish here, we require 5''(0) = 0. In fact, wc expect the free energy to be an 
even function of h, with a cusp at /?. = and t < due to the discontinuity of M on this 
line. One more differentiation of the energy will give the specific heat C ~ t'^/vt-'^ = ^-^ 
We also find that the susceptibility is 

X-T^f^ + Z^V^'"^"^^" ■ (13) 
Hq \ yt ) 



For t < and h small and positive. 



Ho 

E ^ Zp^^_tf/yt-i (14) 
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For t = and h > 0, 



M = g(7c/2)h'^/y'^-^ 

HoVh 

E=^g\Tr/2)h'^''-y^^'y- . (15) 
dTo 

From these expressions we get the following intuition about the scaling free energy g{6): 

1. g{0) controls the singular part of the energy for T > T^. 

2. g{TT) controls the singular part of the energy for T < Tc. 

3. g'{7i/2) controls the energy for t = 0, h ^ 0. 

4. linie^^ g'{9) = — g'{9) controls the expectation value of M for T < T^.. 

5. 5'(7r/2) controls M for t = and /i 7^ 0. 

Here it is convenient to choose Tq and Hq so that HoM(t = 0,h) = h^^^ and HQM(t < 0,h = 
0) = (— t)^. The normalization conditions on t and h then require that g{7r/2) = Uh/d and 
g'ij^) = —1. When it is necessary to distinguish, we will call Hq for the "/(x)" and ''g{9)" 
forms Hf and Hg respectively. Similarly we distinguish Tf and Tg. They are related by 
Hf = Hl^^ and Tf = TgHy^. 

Simulations 

Monte Carlo simulations were run on 16'^, 24^, 32^, 40^, 48^ and 64^ lattices using a 
multiple cluster updating algorithm |jl4 . 

To use a cluster updating algorithm with a nonzero magnetic field, just imagine that in 
addition to the regular bonds with strength J connecting neighboring spins, each spin is 
connected to a fake "magnetizing spin" by a bond of strength H, as illustrated in Fig. ^ 
Then break both "J bonds" and "H bonds" and update clusters according to the usual cluster 
algorithm^]. The "magnetizing spin" is a member of a cluster, and is reflected just like 
any other spin. When evaluating the magnetization of the lattice, we take the components 
of the lattice spins in the current direction of the magnetizing spin. 

Results for the magnetization of the 0(4) model are plotted in Fig. ^. Results from the 
largest lattice size run at each point are shown. The remaining finite size effects are about 
the same size as the statistical error bars. Then, in Fig. ^ the results for h = 0.002, 0.005, 
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0.01 and 0.02 are plotted in the form in Eq. |^. Here we used the values for the critical 
coupling and exponents from Kanaya and Kaya0. 



Then we fit the magnetization results to find an approximate scaling function g{6) for 
the free energy in the form in Eq. |^. A simple parameterization of g{6) which satisfies the 
normalization conditions on g{7c/2) and g'{TT) is 



g{e) = W3 + 2(cos(e/2)-/l72) 
+ao cos(6') 



+ai(cos(3^/2) + 3cos(0/2) -2^1/2) 
+a2{cos{29) + 1) 



+a3(cos(5^/2) - 5 cos(0/2) + 6^1/2) 

(16) 



(c? = 3 in this equation.) 



In this fit I used the energy and magnetization for 0.89 < J < 0.99 and H = 0.005 and 
0.002. The free energy also included an analytic part /a = CH2h^ + Cjit + Cj2t^ + Cj^t'^. 
The resulting g{9) is plotted in Fig. ||. Tg and Hg were 0.44 and 1.31 respectively. The 
magnetization corresponding to this free energy is also plotted in Fig. |^. In principle, the 
critical exponents yt and yh, and the critical coupling are also parameters in this fit. However, 
to get these parameters to the same accuracy as has already been done by Kanaya and 
Kaya[^ or Butera and Comi|]ll| would require a careful correction for finite size effects. 



and care in using only data for small enough t and h that corrections to scaling are small. 
Therefore, the exponents and critical coupling were fixed to those found by Kanaya and 
Kaya. Because of the remaining finite size effects and corrections to scaling, the of this fit 
was very bad (243 for 30 degrees of freedom). However, since the results are already many 
times more accurate than the QCD data with which we intend to compare, there is little 
incentive to make the necessary corrections. This scaling function was then converted to 
the "/(x)" form (by computing the resulting magnetization as a function of t for h = 0.002 
and plotting according to Eq. |^.), and plotted as a solid line in Fig. |, where it can be seen 
to describe the magnetization quite well. (It is necessary to convert the normalization of 
if: and h used in the "(y'(^)" form to the conventional normalizations for the "/(x)" form: 
Hf = H^+^ and Tf = TgHy^.) In this fi gure I have also included the mean field form of the 
scaling function and the epsilon expansion form. 

Fig. ^ also shows the properly normalized asymptotic form for f{x) clS X ^ CX^ 5 X 
It can be seen that the scaling function approaches this asymptotic form quite slowly. This 
is the region where the long distance physics is dominated by the Goldstone bosons. In 
particular, we expect that for t < and h small, the magnetization takes the form M = 
M{t,0) + A/i^/^, so that the susceptibility diverges at /i = for all t < 0[T^. The fitting 
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Figure 2: Cluster updating with a magnetic field. A single "magnetizing spin" is coupled 
to every spin in the lattice by a bond with strength H. 

function Eq. ^ should really be modified to support this behavior at = tt, but this problem 
seems to occur in a region beyond where these 0(4) results, and the QCD results to which 
they will be compared, are taken. 

The Monte Carlo scaling function and the epsilon expansion are in good agreement for 
t < because of the normalization condition on t, M{t < 0, /i = 0) = (— t)^. Had we chosen 
the equally sensible normalization condition X = ^ = for t > and = 0, we would 
have found agreement of the Monte Carlo and epsilon expansion for t > with a discrepancy 
for t < 0. 

Acknowledgement s 

Tom Blum, Carleton DeTar and Bob Sugar have contributed greatly to this work. I thank 
Joe Rudnick for a valuable conversation, particularly for pointing out to me the diverging 
susceptibility along the —t axis. This work was supported by DOE grant DE-FG03-95ER- 
40906. 



References 

[1] R.D. Pisarski and F. Wilczek, Phys. Rev. D 46, 4657 (1992); F. Wilczek, J. Mod. 
Phys. A7, 3911 (1992). 



7 



0.5 ^ \ I I I ^ \ \ I I I I \ \ I I 

I 0(4) magnetization ^ ^'i 

0.4 — Q □ - o ^ 

o ° □ ^'^ S * 

0.3 — Q ° ;k ^'^ o o * 



Q Ho 

'G □ 



o o ^ 
o ^ 



0.2 — H ^ ^ 

□ ^ 



□ 



□ □ ^ 



0.1 K Q ° - > 

G 

_G ^ 



0.0 



J I I \ I I I I I L 



0.85 0.90 0.95 1.00 

J 
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Figure 4: 0(4) magnetization for h = 0.002, 0.005 and 0.01 plotted as a scaling function 
in the conventional form. I use the results of Kanaya and Kaya0 for the critical coupling 
and critical exponents. In this plot the points for h = 0.05 are plotted with plusses, those 
for 0.02 with crosses, 0.01 with diamonds, 0.005 with octagons and 0.002 with squares. Also 



shown are the asymptotic forms f{x) 



'— x)^ as X 



—oo (a) and f{x) ~ Cx as s — > cxd 



(b). A four parameter fit to the scaling function from the ^^g{0)" form is shown with a solid 
line (c), running through the entire graph. The mean field scaling function, and the second 
order epsilon expansion scaling function [|T^ are shown, labelled "e" and "mf " . 
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Figure 5: Four parameter fit to the scaling function g{9) for the 0(4) free energy. 
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Figure 6: The scaling part of the magnetization for 0(4), cor- 
responding to the scahng function in Fig. ^. 
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